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ADAPTIVE  METHODS  AND  ERROR  ESTIMATION  FOR 
ELLIPTIC  PROBLEMS  OF  STRUCTURAL  MECHANICS 

I.  Babuska,*^  A.  Miller*  and  M.  Vogelius* 


Abstract.  We  develop  a  class  of  feedback  finite  element  methods 
for  a  one  dimensional  boundary  value  problem,  and  we  state  some  mathe¬ 
matical  results  about  their  optimality  with  respect  to  various  perfor¬ 
mance  measures.  Analogous  ideas  and  results  are  discussed  for  a  two 
dimensional  problem  from  plane  elasticity.  A  post-processing  technique 
aimed  at  computing  the  value  of  a  functional  with  high  accuracy  is  also 
described  and  illustrated  by  an  example. 

1.  Introduction.  The  finite  element  method  is  today  the  main  tool 
in  computational  structural  mechanics,  and  the  reliability  of  its  re¬ 
sults  is  thus  of  major  importance.  Judging  the  reliability  of  an  en¬ 
gineering  computation  is  a  complex  matter  involving  many  considerations. 
However,  we  shall  only  concern  ourselves  here  with  the  accuracy  of 
the  finite  element  solution  Itself,  that  is,  we  compare  it  with  the 
exact  solution  of  the  chosen  mathematical  model.  The  goal  is  typically 
either  to  obtain  an  approximation  with  an  accuracy  of,  say  5-10% 
measured  in  a  suitable  norm,  or  to  find  values  of  the  displacements, . 
stresses,  stress-intensity  factors  at  specific  points,  etc.,  with,  say,  1% 
accuracy. 

One  of  the  main  decisions  in  a  finite  element  computation  is  the 
choice  of  the  mesh.  In  practice  various  mesh  generators  are  used; 
however,  a  proper  choice  of  mesh  should  be  closely  related  to  the  aim 
of  the  computation.  An  effective  mesh  cannot  be  constructed  without 
some  feedback,  either  in  a  user  Interactive  fashion,  or  preferably 
through  a  largely  automated  procedure. 

In  this  paper  we  discuss  the  feedback  finite  element  method  and 
adaptivity  in  relation  to  given  performance  measures.  We  address  these 
Issues  for  two  different  computational  goals: 

a)  to  obtain  an  approximate  solution  with  prescribed  accuracy  in 
the  energy  norm, 

b)  to  obtain  an  approximate  value  of  a  functional  (the  stress  value 
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at  a  particular  point)  with  a  prescribed  error-tolerance. 

In  Section  2  we  discuss  for  illustration  a  simple  one  dimensional 
model  problem.  We  introduce  the  notion  of  a  feedback  finite  element 
method  and  define  adaptivity  to  mean  optimality  with  respect  to  a  given 
performance  measure.  The  Theorems  2.2,  2.3  and  2.4  verify  adaptivity 
for  some  specific  feedback  finite  element  methods. 

Section  3  contains  a  discussion  of  the  adaptive  features  of  the 
FEARS  program  (Finite  Element  Adaptive  Research  Solver)  developed  at 
the  University  of  Maryland.  The  numerical  example  is  the  computation 
of  the  stress  state  of  a  cracked  panel. 

In  Section  4  we  address  the  question  of  adaptivity  in  a  case  where 
the  aim  of  the  computation  is  to  find  the  maximal  stress  of  a  twisted 
bar. 

2.  The  One  Dimensional  Problem.  Consider  the  model  two  point  bound¬ 
ary  value  problem  on  I  ■  (0,1) 

(2.1a)  Lu  5  (a(x)  ^jj)  +  b(x)u  -  f(x),  x(  I 

(2.1b)  u(0)  -  u(l)  -  0 

where  a,b  €  1^(1)  with 

(2.2a)  0  <  aj  <  a(x)  <  cu  <  “ 

(2.2b)  0  <  b(x)  <  <*2  <  00 

and 

(2.2c)  f  €  H_1(I) 

°1 

The  solution  u  i  H  (I)  is  to  be  understood  in  the  usual  weak  form 

B(u,v)  -  f  fvdx  V  v  €  H1(I) 

J0 

where  B  denotes  the  bilinear  form 

(2.3)  B(u,v)  «  (a  ^  ~  +  buv)dx  u,v  (  H*(I). 

q  qx  ax 

The  finite  element  solution  is  Introduced  in  the  standard  way.  Let 
A  be  an  arbitrary  mesh 
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A  A  A 

°  -  x0  "  X1  < 


and  let  1^  •  *  x^  -  x^  ,  j  -  1,...,N(A),  h(A) 


max  h. 


The  Xj  will  be  called  nodal  points,  the  Ij  elements.  S^(A) 

°1 

denotes  the  subspace  of  H  (I)  consisting  of  functions  that  are  poly- 

A 

nomials  of  degree  5  p  on  each  element  Ij .  The  finite  element  solu¬ 


tion  u(A)  €  SP(A)  is  defined  by 


(2.4) 


B(u(A) ,v)  -  B(u 


.v,  -  I1 


Vv  €  SP(A) , 


i.e.,  u(A)  is  the  elliptic  projection  of  u  onto  SP(A).  The  expres¬ 
sion  e(A)  «  u  -  u(A)  is  the  error  of  the  finite  element  solution. 
Define  IMIg  to  be  the  energy  norm  of  v: 

livllg  -  [B(v,v)]1/2  . 

We  are  interested  in  effectively  computing  u(A)  so  that  ||e(A)||^, 

or  the  relative  error  |je(A)||-,/||u|L,  is  smaller  than  some  given  tol- 

£  £ 

erance. 

For  any  subinterval  I  $  I,  let  R(I-I)  denote  the  elliptic  pro¬ 
jection  of  onto  the  space 

(v  (  H^I)  :  v (x)  -  0  V  x  €  1-1} 

and  define  the  error  indicatiors  n^dj)  *  Hj,  j  *  1,...,N(A)  and 
the  error  estimator  e(A)  as  follows: 


(2.5) 


(2.6) 


lR(I-Ij)e(A)|lE 


p(A)  .  jV 

A  ** 


The  error  indicators  can  be  estimated  by  expressions  involving  the 
residuals  r,  *  (f  -  Lu(A))|  ,,  see  e.g.  (8). 

j  I  T/l 

J 
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We  underline  the  fact  that  the  error  indicators  and  the  estimators  are 
directly  related  to  the  goal  of  the  computation,  in  this  case  to  achieve 
the  desired  accuracy  measured  in  the  energy  norm. 

An  important  result  is 

Theorem  2.1.  There  exist  C  <  00  depending  on  o^f  a2  of 
(2.2  a,b)  but  independent  of  a,  b,  f  and  A  such  that 

(2.7)  e(A)  S  ||e(A)||E<  Ce(A). 

If  we  make  additional  regularity  assumptions  about  a,  b,  f  then 
it  is  possible  to  prove  that  the  effectivity  index 

(2.8)  6(A)  -  )|e(A)f|E  "  1  as  He(A)liE  ■*  ”  » 


for  more  details,  see  [8]. 

Let  us  remark  that  e(A)  is  not  in  general  an  upper  bound  for  the 
error  (because  C  >  1).  It  is  possible  to  estimate  C  in  (2.7),  but 
the  estimate  we  have  is  very  pessimistic.  It  is  also  possible  to  find 
a  computationally  more  expensive  estimator,  that  provides  an  upper 
bound  (cf.  [8]). 


Let  us  now  define  a  feedback  finite  element  method.  In  very  general 
terms  it  consists  of  consecutive  refinements  of  the  meshes  in  dependence 
on  computed  results.  More  precisely  in  our  model  problem  we  consider 
only  binary  meshes,  i.e.,  the  nodal  points  are  of  the  form  J2—  ,  and 
the  transition  operator  performs  a  bisection  of  certain  elements.  By 
successive  application  of  the  transition  operator  a  sequence  of  meshes 
is  created  and  the  corresponding  finite  element  solutions  are  computed, 
until  the  process  is  terminated  through  a  stopping  criterion. 


The  sequence  of  binary  meshes  A^.Aj,.*.  (with  S^(A^)  £  S**(A^+j) 
together  with  the  u(A^)  is  called  £  trajectory.  A  feedback  finite 

element  method  creates  a  trajectory  which  depends  on  the  exact  solution 
u,  the  coefficients  a,  b  and  the  initial  mesh  A^.  Given  a  perfor¬ 
mance  measure  on  the  set  of  trajectories  a  feedback  finite  element  me¬ 
thod  will  be  called  optimal  or  adaptive  relative  to  this  measure  and  the 
set  H  if  for  any  (u,a,b,A^)  (  H  it  creates  a  trajectory  with  optimal 

performance  measure  value  (see  [11]). 


Below  we  define  three  simple  transition  operators 
r  3 1 

A1  ;  in  each  case  we  list  the  elements 
*  Ai+1! 


*m. 


m 


sing  from  A.,  to  A., 


that  are  bisected  in  pas- 
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m  a4  t 

l  :  all  elements  I.  for  which  n.  >  Bmaxn  (I) 
i  J  J  I€A, 


0  5  B  5  1 


12]  A 


l  !  I j (!)*’“ *Ij (k) *  where  nj(l)  “  nj( 


(2)  nj(N(A1)) 


in  an  ordering  (by  magnitude)  of  the  indicators  and 


k  «  tvN(Ai)  ]  +  1 


0  5  y  <  1 


([a]  denotes  the  integer  part  of  a). 

13]  [1]  [2] 

l  >1  u  l  • 

ill 

The  transition  operator  A  ^  bisects  all  elements  whose  indicators 
are  larger  than  or  equal  to  6  times  the  largest  indicator;  for  6  «  0 


this  implies  bisection  of  all  elements.  A 


bisects  approximately 


the  fraction  y  of  the  elements  and  the  transition  operator  A1  J  bi¬ 
sects  all  elements  with  indicators  larger  than  or  equal  to  B  times 
the  largest  one,  and  at  least  the  fraction  y. 

In  practice  it  is  advantageous  to  use  a  more  composite  transition 

A 

operator.  Define  the  predictor  n(I^)  by 


fi(lj  )  -  min 


Ai  Ai  2  \ 

In  il,1) l  A.  a .  \ 

7$r- "  "4 


0  <  x  <  1 


where  I?  (  A.  ,  k  <  i,  is  the  direct  predecessor  of  I.  .  The  direct 

Ai  Ak  j 

predecessor  of  I  is  defined  as  the  I~  (with  maximal  k)  that 

A.  A.  Aj  ^  A. 

satisfies  Ij  ?  1^  ,  h^  ■  2h  ,  except  for  I  (  A ^  where  the 


direct  predecessor  is  defined  to  be  I 


?i  is  called  a  predictor 


since  it  uses  past  experience  to  predict  the  effect  of  bisection  on 

a  tral  na  a  f  f  Ua  4  4  a  aw  TU  a  f  «•  . .  i  [ A]  u . . 


the  value  of  the  Indicator.  The  transition  operator  A 
Ai 

all  elements  I  that  satisfy 


bisects 


5 
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i(I.i)  >  6  max  fi(I) 
3 


0  5  6  5  1. 


The  transition  operator,  A 


[5] 


[4] 


used  in  the  FEARS  program  is  a  di¬ 


rect  extension  of  the  operator  A  ,  including  the  so-called  short 
passes:  whenever  the  increase  in  the  number  of  elements  is  less  than  a 
preset  percentage  of  the  total  number  of  element  of  A^  no  finite 

element  solution  is  computed  on  A..,.  Instead  indicators  for  the  ele- 

i+1  a 

ments  of  A^+^  are  constructed  directly  from  the  fi(Ij  )  and  the 


operator  A 


[4] 


is  applied  to  A 


1* 


• A^ ,A^+^  with  this  set  of  indica- 


i+1  ui+l 

n  (Ij  ) .  The  process  continues  until  the  increase  in  the 


tors  for 

number  of  elements  is  larger  than  the  preset  percentage  of  the  total 
number  of  elements  of  A^,  at  which  point  the  short  pass  is  concluded 
and  the  finite  element  solution  is  computed  on  the  resulting  mesh 

The  transition  operators  A^  all  have  the  property  that 

A.  A  j 

they  bisect  at  least  one  element  with  n.1  >  B  max  n  x(l)»  0  <  6  <  1 

J  ~  I€A^ 

(Indeed,  they  all  bisect  the  one  with  the  largest  indicator).  We  shall 
call  a  transition  operator,  and  the  corresponding  feedback  finite  ele¬ 
ment  method,  with  this  property  regular .  The  following  theorem  is 
proven  in  [ 9 ] , 


Theorem  2.2.  Let  u  €  H  (I)  denote  the  solution  of  (2.1  a,b)  and 
let  u(A^)  i  *  1,2,...  be  the  sequence  of  solutions  created  by  a 

regular  feedback  finite  element  method,  then 


I  •<a1)||e 


0  as  i  -►  », 


We  note  that  h(A^)  does  not  necessarily  converge  to  zero. 

If  we  define  the  performance  measure  (of  a  trajectory)  to  be  1 
whenever  u(a^)  converges  to  u  (in  the  energy  norm)  and  0  other¬ 
wise,  then  Theorem  2.1  simply  states  that  any  regular  feedback  finite 
element  method  is  adaptive  with  respect  to  this  measure  and  the  set  of 

°  j 

u  £  H  (I),  a,  b  satisfying  (2.2  a)-(2.2  c)  and  arbitrary 
Let  $  be  defined  as  follows 


$(u,N)  *  inf  Inf  ||u  -  vjl 

.' :  binary  mesh  v^S^(A) 

n(a)*n 


Adaptive  Methods  and  Error  Estimation 


and  consider  the  performance  measure  that  is  1  whenever 
(2.9)  || u  -  u(Ai)|iES  CWu.NCAj)) 

(for  some  C  Independent  of  i),  and  0  otherwise. 

In  [ 9  ]  we  analyze  the  adaptivity  with  respect  to  this  measure  in 
the  case  when  a  in  addition  to  satisfy  (2.2  a)  is  Holder  continuous 

with  exponent  *1/2.  The  transition  operators  A  ^  -A  ^  are  not 
adaptive  with  respect  to  this  measure  for  all  u  €  H*(I).  However, 
it  is  possible  to  verify  adaptivity  for  a  rather  large  class  of  func- 

O  1 

tion  K  9  H  (I).  We  shall  not  give  here  the  exact  definition  of  K, 
instead  we  refer  to  [  9  ] .  The  class  K  contains  for  instance  all 
functions  that  are  smooth  or  have  singular  behaviour  as  xa,  a  >  1/2 

(xa-x  (  K) .  The  following  theorem  is  proven  in  [9  ]: 

Theorem  2.3.  Assume  that  a  is  Hblder  continuous  with  exponent 
1/2  and  a,  b  satisfy  (2.2  a,b).  Let  A^  be  a  sequence  of  meshes 

generated  by  a  feedback  finite  element  method  based  on  either  of  the 

transition  operators  A^  or  A^  with  0  <  6;  then 

||e(A1)||E<  C(u)4»(u,N(A1)) 


provided  u  (  K. 

Theorem  2.4.  Let  the  assumption  be  as  in  the  previous  theorem.  Let 

A.  be  a  sequence  of  meshes  generated  by  a  feedback  finite  element 

1  r  2 1  r  3 1 

method  based  on  either  of  the  transition  operators  A  or  A1  J 

with  6  >  0. 

For  any  u  €  K  there  exists  Yq  such  that  if  0  <  y  <  Yq»  then 

||e(Ai)||Ei  C(u)<t>(u,N(Aj))  • 

It  should  be  noted  that  A^  and  A ^  with  8*0  do  not  in 

general  (e.g.  when  u  *  xa  -  x)  produce  optimally  convergent  finite 

element  solutions;  consequently  A^  and  A^,  with  8*0,  are 
not  adaptive  relative  to  K  (they  are  adaptive  relative  to  a  sub¬ 
class  of  smooth  functions  c  K) . 

We  also  observe  that  in  most  practical  cases  CjN  ^  5  <t>(u,N)<C2N 

whether  u  is  smooth  or  not  (e.g.  for  u  *  xa-x,  a  >  1/2).  Thus 
a  feedback  finite  element  method  which  is  adaptive  with  respect  to 
this  performance  measure  leads  to  a  rate  of  convergence  that  (in 


9  •  •  •  4 


V  ^ 


If 
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practical  cases)  is  not  influenced  by  the  smoothness  of  the  solution. 

Let  the  cost  of  computing  the  finite  element  solution  u(A^),  the 
Indicators  rij*»  the  estimator  e^)  and  creating  the  mesh  A^+1  be 
denoted  by  (p(A^).  The  total  cost  ^(A^)  =  <p(A^)  and  the 

proportionality  factor 


e  (At) 


<('(Ai) 

cpCAj^) 


describes  the  effectivity  of  the  feedback  method.  Assuming  that 

\  r 2 ] 

<p(A^)  *  N  (A^)  ,  X  >  1,  then  for  the  case  of  the  operators  A 

A^  we  have 


and  hence 


B(V  5  1+7  "<4i+i> 


j-i  .  1 1 

nA.)  5  l  ( - i-j  )  NA(A.) 

3  i-o\(i+y)a/  3 


,_0±lL_NX(  j. 

(l+y)*-l  3 

r  2 1  r  31 

Thus  for  the  operators  A  and  A  we  get 

tu.»  5  . 

ci+f)  -i 

A  similar  inequality  is  not  (in  general)  true  for  A^  or  A^; 
whenever  only  very  few  new  elements  are  added  in  each  transition  the 

ratio  C(A^)  becomes  unbounded  as  i  -*•  °°.  The  operator  A^  is 

created  exactly  in  order  to  avoid  this  difficulty. 

In  the  preceding  analysis  we  considered  only  the  case  of  one  given 
load  function  f,  but  the  generalization  to  multiple  loads  is  straight¬ 
forward:  We  Interpret  the  multiple  load  problem  as  a  system  of  equa¬ 
tions  for  u  ,  1  i  l  f  n.  We  introduce  the  energy  norm 
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n  ,11/2 

K  -  W'eJ 

and  we  define  the  indicators  and  estimators  by  the  appropriate  projec¬ 
tions  in  this  norm. 

3.  Two  Dimensional  Problems:  The  FEARS  Program.  Many  of  the  ideas 
discussed  for  the  simple  one  dimensional  problem  of  the  previous  section 
have  been  implemented  in  a  two  dimensional  setting  in  the  FEARS  program. 
FEARS  is  able  to  handle  symmetric  elliptic  systems  of  two  equations  on 
regions  which  are  unions  of  curvilinear  quadrilaterals.  Within  the 
program,  each  of  these  quadrilaterals  is  mapped  onto  a  unit  square.  The 
actual  finite  element  calculations  are  carried  out  on  these  mapped 
squares  using  conforming  square  bilinear  elements.  A  typical  FEARS 
mesh  for  the  shaded  region  of  Fig.  1  is  shown  in  Fig.  2.  In  this  paper 
we  shall  illustrate  some  of  the  major  features  of  the  program  by  refer¬ 
ence  to  a  model  problem  and  by  analogy  with  the  one  dimensional  case 
discussed  in  Section  2.  For  further  details  and  a  fuller  discussion 
see  [1],  [5],  [7]  and  [10]. 

As  the  model  problem,  consider  the  stress  state  of  the  cracked  panel 
illustrated  in  Figure  1.  We  suppose  that  plane  strain  assumptions 
apply  with  v  «  Poisson’s  ratio  *  .3  and  E  *  Young's  modulus  •  1. 
(Because  of  the  symmetries  which  are  present,  we  need  only  compute  on 
the  shaded  quarter  panel.)  As  in  Section  2,  we  shall  be  interested  in 
judging  the  accuracy  of  our  finite  element  approximation  in  terms  of  the 
strain  energy  norm  of  the  error.  The  exact  solution  of  this  model  pro- 
1/2 

blem  has  an  r  -type  singularity  in  displacements  near  the  crack  tip, 
as  well  less  severe  singular  behaviour  at  the  corners.  Indeed,  the 

solution  is  in  the  Besov  space  B^  for  u  *  3/2  (and  no  greater) . 

-1/4 

As  a  consequence,  using  uniform  meshes  we  would  only  expect  an  0(N  ) 

rate  of  convergence  in  the  energy  norm.  Better  convergence  rates  than 
this  can  be  achieved  by  appropriate  mesh  refinements. 

As  in  the  one  dimensional  case,  error  Indicators  can  be  associated 
with  each  element  of  the  mesh.  However,  the  computation  of  such  an 
indicator  now  needs  information  not  only  about  the  finite  element  solu¬ 
tion  on  the  element  itself,  but  also  on  adjacent  elements.  Error  esti¬ 
mators  are  constructed  from  these  Indicators  exactly  as  in  Section  2. 

The  transition  operators  used  are  direct  analogues  of  A^  -  A^. 
Theorems  2.1  and  2.2  carry  over  directly  to  the  present  case,  while 
(2.8)  holds  under  some  additional  assumptions  on  the  meshes.  Experi¬ 
ments  seem  to  indicate  that  these  extra  assumptions  are  necessary  in 
the  two  dimensional  case,  and  that  the  analogues  of  Theorems  2.3  and  2.4 
appear  to  hold. 

In  Figure  3  we  have  plotted  the  relative  error  in  the  energy  norm 
against  the  number  of  degrees-of-freedom  N  for  two  feedback  trajec¬ 
tories  produced  by  FEARS:  I  was  created  using  the  transition  operator 
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with  P  *  1  and  y  =  .2;  II  used  A^  with  6  =  0  (which 
leads  to  a  uniform  mesh).  As  theory  would  predict  the  trajectory  II 

-1/4 

of  uniform  meshes  shows  an  0(N  )  rate  of  convergence.  However  I 


-1/2 

shows  an  0(N  )  rate.  This  is  the  maximum  possible  rate  for  ele¬ 

ments  of  degree  1.  The  mesh  shown  in  Fig.  2  corresponds  to  the  last 
computed  mesh  of  I#  and,  as  expected,  it  exhibits  a  severe  refinement 
around  the  crack  tip.  The  ratio  of  maximum  element  size  to  minimum 
element  size  for  this  mesh  is  256. 


Rather  more  meaningful  than  the  dependence  of  relative  error  on  N 
is  its  relation  to  total  cost  (i.e.  computer  time).  Fig.  4  is  a  plot 
of  accuracy  against  computer  time  for  the  trajectories  I  and  II,  and 

T4 1 

for  the  trajectory  III  which  was  created  by  A  with  6=1 
i.e.  no  Mshort  passes11  were  permitted.  Again  we  see  the 
marked  superiority  of  I  over  II.  Also,  for  a  given  accurary,  I  re¬ 
quires  approximately  half  as  much  time  as  III.  The  trajectory  III 
"wastes”  time  solving  on  meshes  that  only  differ  slightly  from  one  ano¬ 
ther.  The  proportionality  factors  for  trajectory  I  are  in  the  range 
1.5  <  <  2.5,  while  those  for  III  lie  in  the  range  3  <  £  <  5.  For 

problems  with  more  severe  singular  behaviour  than  our  model  problem 
(e.g,  a  clamped-free  crack)  the  relative  economy  of  I  over  III  could  be 
expected  to  be  even  greater. 

Finally  in  Figure  5  we  show  the  dependence  of  the  effectivity  index 
6  on  solution  accuracy.  The  plot  is  consistent  with  6  1  as 

|| e || ^  -►  0.  Let  us  just  mention  that  FEARS  has  two  kinds  of  error  indi¬ 
cators  and  estimators.  Asymptotically  both  behave  similarly,  though 
one  is  always  greater  than  the  other.  In  our  calculations  we  used  the 
larger.  This  is  usually  preferable  since  it  provides  extra  "safety" 
in  the  stopping  criterion. 
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Figure  5 


4.  Post-Processing  Calculations.  So  far  we  have  only  been  judging 
the  accuracy  of  a  finite  element  solution  in  terms  of  the  energy  norm 
|| u  -  u(A)|L  of  the  finite  element  error*  However,  in  practice,  often 

£j 

the  goal  of  a  finite  element  calculation  is  not  to  obtain  high  accuracy 
in  the  energy  norm  as  such,  but  rather  to  find  sufficiently  accurate 
approximations  for  values  of  certain  functionals  (for  instance,  stress 
at  a  point,  stress  intensity  factor  at  a  crack  tip,  average  flow  rate 
through  a  surface,  etc.).  In  this  section  we  shall  address  this  point 
by  reference  to  a  model  problem.  For  a  more  complete  discussion  see 
[2],  [3]  and  [4]. 

Let  us  consider  a  stress  function  formulation  of  the  torsion  problem 
for  a  square  bar.  The  governing  boundary  value  problem  is 

(4.1a)  V^u  «  -1  on  ft  *  (-1,1)  x  C— 1 » 1) 

(4.1b)  u  =  0  on  the  boundary  9ft  of  ft. 

We  shall  be  interested  in  this  problem  only  as  a  means  of  finding  the 
(maximum)  stress 


(4.2) 


au 

ax. 


x*(l,0) 


Given  the  finite  element  approximation  u(A),  the  most  common  approach 

would  be  to  use  \.P(A)  ag  an  approximation  to  a •  However, 

X1  x*(l,0) 

we  shall  employ  a  post-processing  technique  which  yields  better  results 
than  this  standard  polntwise  evaluation  approach.  It  may  be  shown  that 
for  the  problem  (4.1) 


(4.3) 


_au_ 

ax. 


x*(l ,0) 


(u^  +  e  )  d  A , 


where 
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(4.4b) 


. .  i  *lL  _ , 

’  <*,-!>  o 


and  (pQ  is  any  smooth  function  which  ensures  that  a  0  on 
3ft  -  (1,0).  For  instance,  we  could  take 


(4,4c) 


(x^ir+i 


X  -1  X--1 

— —  +  — — 
5  2x/ 

x2+4 


The  relation  (4.3)  can  readily  be  verified  by  an  integration  by  parts 
and  a  simple  limiting  argument  (see  [2]).  Obviously  (4.3)  suggest  that 


(4.5) 


(u(a) 
J  o 


(A)£  +  <p)dA 


as  an  approximation  for  o.  The  expression  (4.5)  is  one  example  of  an 
extraction  expression  for  o,  and  is  the  basis  of  the  post-processing 
calculation  that  we  shall  use.  Let  us  just  mention  that  the  standard 
pointwise  evaluation  approach  can  also  be  written  if  the  form  (4.5)  if 

1  Xl~^ 

we  formally  let  cpn  *  — - = — tz  .  (So  cp  =  0  and  £  is  a  deriva- 

0  (xrl)Z+x‘ 

tive  of  the  Dirac  delta  function.) 

Let  us  now  discuss  the  accuracy  of  the  post-processed  value  a (A). 
Clearly 


(4.6) 


o  -  o  (A) 


* 

(u  -  u(A)KdA, 


and  if  we  now  introduce  an  auxiliary  problem 


(4.7) 


V >  -  -f. 


on  n 


on  9ft 


then  (4.6)  gives 


(4.8)  o  -  a  (A)  *  V(u  -  u(A))  •  V (0  -  4/(A))dA. 

Jft 

This  shows  that  the  accuracy  of  a (A)  depends  not  only  on  the  accuracy 
of  u(  )  but  also  on  the  error  (\i»-^( A)).  Clearly  the  choice  of 
finite  element  mesh  should  recognize  this  fact. 

There  are  some  alternate  forms  of  (4.8): 
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(4.9)  o  -  a  (A)  *  ^(!|(u  + iJj)  -  (u  +  (A)  ||g  -  ||(u-ij/)  -  (u-4>)(A)!|^) 

and  somewhat  less  sharply 

(4.10a)  |o  -  o(A)  |  5  ||u-u(A)||e  ||  -  (|»(A)He 

(4.10b)  |o-o(A)|  -  ^  (||u-u(A)||2  +  a2l!^-iHA)||^)> 

where  for  any  function  v  we  have  used  ji v ||  =  (i  |vv|  2dA)1/2  to 

b  Jft 

denote  the  energy  norm  associated  with  (4.1).  Let  us  mention  that 
(4.10a)  and  (4.10b)  in  general  give  upper  bounds  for  the  error  in  a (A) 
as  they  fail  to  take  account  of  any  cancellation  in  the  integral  (4.8). 

The  feedback  finite  element  method  already  discussed  in  Sections  2 
and  3  can  be  applied  to  this  post-processing  calculation  provided  the 
appropriate  performance  measures  are  chosen.  It  is  possible  to  prove 
that  Theorem  2.2  (relating  to  the  convergence  measure)  holds  if  the 
feedback  method  is  regular  with  respect  to  the  indicators  associated 

2  22  1/2 

with  the  norm  (||uj|E  +  a  IMIg)  ]  .  (See  the  remarks  on  multiple  load 

computations  at  the  end  of  Section  2.)  However,  it  seems  too  much  to 
ask  that  the  direct  analogues  of  Theorem  2.3  and  2.4  should  hold. 
Because  of  the  possibility  of  cancellation  in  the  integral  (4.8)  we  can 
expect  that 


4>(N)  -  inf  |o  -  o(A)  | 

A 

N(a)«N 

could  be  quite  small  (possibly  zero),  even  for  small  values  of  N. 
However,  we  can  frame  a  less  demanding  performance  measure  based  on 
(4.10b).  Let 

*<N)  -  inf  || u  -  u(A)||2  +  c2||*  -  *(A)||2  . 

A  E  E 

N(A)«N 

The  analogues  of  Theorems  2.3  and  2.4  are  now  closely  related  to  the 
corresponding  results  in  Section  3. 

A  posteriori  estimates  for  this  post-processing  calculation  may  be 
based  on  (4.9)  and  (4.10).  FEARS  gives  estimators  c^(A),  Cj( A)  for 

||(u  +  ifi)  -  (u  +  U;)(A)||E  and  ||(u-ty)  -  (u  -  )  (A )  H  E  respectively  with 

effectivity  indices  which  converge  to  1.  Now  taking  the  obvious  arith¬ 
metic  combinations  of  there  estimators  suggested  by  (4.9), 
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(4.11)  e  (A)  -  •|(c1(A)2- e2(A)2), 

gives  an  estimator  for  the  error  o-o  (A)  in  the  post-processed  value* 
This  estimator  will  have  an  effectivity  index  which  converges  to  1, 
provided  the  subtraction  in  (4*11)  does  not  lead  to  significant  cancel¬ 
lation.  If  this  critical  case  occurs,  we  cannot  expect  e(A)  to  be 
reliable.  We  may  nonetheless  use  estimators  based  on  (4.10).  Such 
estimators  however  will  be  pessimistic,  especially  in  this  critical 
case.  See  [4]  for  further  discussion  of  this  point. 

Table  1  shows  some  numerical  results  for  our  model  problem.  The 
problems  (4.1)  and  (4.7)  were  run  as  a  multiple  load  system  using  FEARS 
Since  both  solutions  u  and  \J/  are  rather  smooth  it  is  not  surprising 
that,  at  least  initially,  all  the  meshes  produced  were  uniform.  Notice 
the  following  features  of  the  results  in  Table  1: 

(a)  |E u  -  u(A)  ||  displays  the  usual  0(h)  rate  of  convergence 

£j 

expected  of  bilinear  elements  for  a  smooth  exact  solution  u. 

(b)  The  effectivity  index  of  the  estimator  for  ||u  -  u (A)  ||  is 

close  to  1. 

(c)  The  post-processed  value  a (A)  is  markedly  superior  to 

1  .  The  standard  value  is  converging  as  0(h),  whereas  the 

1  >x-(l,0) 

2 

convergence  rate  of  the  post-processed  value  is  0(h  ) . 

(d)  The  effectivity  index  for  the  estimator  (4.11)  is  very  close  to 
1,  whereas  for  the  estimator  based  on  (4.10a)  it  seems  to  stabilize 
around  1.3.  The  fact  that  (4.10a)  leads  to  an  estimator  which  is  only 
a  slight  overestimate  of  |o  -  o(A)|  indicates  that  there  is  little 
cancellation  in  the  integral  (4.8).  The  high  accuracy  in  o(A)  can 

be  explained  solely  by  (4.10a).  Indeed,  for  this  problem,  since  <l>  is 
smooth  we  expect  |lty  -  il>(A)|L  -  0(h).  Thus  (4.10a)  gives 

|o  -  o (A)  ■  0(h  ),  which  is  consistent  with  (c)  above. 
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TABLE  1  i 

( 

Accuracy  of  the  maximum  stresses  calculated  by  post-processing, 
and  a  posteriori  error  estimates. 


NUMBER  OF  ELEMENTS 

(Uniform  mesh,  h  ■  element  size) 

16 

(h* . 5) 

64 

(h® . 25) 

256 

(h*. 125 

ENERGY  NORM  OF  THE  ERROR 
||u  -  u(A)HE/||u|lE 

EFFECTIVITY  INDEX  OF  ESTIMATOR 

FOR  |lu  -  u(A)I|e 

30.1% 

.94 

15.2% 

.96 

7.62% 

.98 

ERROR  IN  STANDARD  VALUE 

(o-8"<4)  >/M 

dXl  x-(l,0) 

29% 

16% 

8.7% 

ERROR  IN  POST-PROCESSED  VALUE 

|o  -  o  (A)  |  / 1 a  | 

1.5% 

.37% 

.089% 

EFFECTIVITY  INDEX  OF 

THE  ESTIMATOR  BASED  ON 

(4.11) 

.99 

1.02 

1.01 

(4.10a) 

1 

1.24 

_ 1 

1.28 

1.29 
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